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Understanding the acid-base behavior of silica surfaces is critical for many nanoscience and bio- 
nano interface applications. Silanol groups (SiOH) on silica surfaces exhibit two acidity constants — 
one as acidic as vinegar — but their structural basis remains controversial. The atomic details of the 
more acidic silanol site govern not just the overall surface charge density at near neutral solution 
pH, but also how ions and bio-molecules interacts with and bind to silica immersed in water. Using 
ab initio molecular dynamics simulations and multiple representative crystalline silica surfaces, we 
determine the deprotonation free energies of silanol groups with different structural motifs. We 
show that previously proposed motifs related to chemical connectivity or inter-silanol hydrogen 
bonds do not yield high acidity. Instead, a plausible candiate for pKa=4.5 silanol groups may be 
found in locally strained or defected regions with sparse silanol coverage. In the process, irreversible 
ring-opening reactions of strained silica trimer rings in contact with liquid water are observed. 
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INTRODUCTION 

Deprotonation of silanol (SiOH) groups [HQ at water- 
silica interfaces is one of the most common and impor- 
tant, yet intriguing, interfacial chemical reactions. Sil- 
ica (Si02) is a major component of rocks and lines the 
channels of many nanofluidic devices. Deprotona- 
tion governs dissolution rates, ^] affects lipid binding to 
silica nanostructures,[3] creates negative surface charges 
that can be tuned vifith moderate changes in solution 
pH to perform desalination and ion gating, Q and may 
even hinder extraction of positively charged crude oil 
components 6] from underground deposits. In particu- 
lar, the atomic level structural details of deprotonated 
SiOH groups govern both the overall surface charge den- 
sity and the binding of ions and molecules to immersed 
silica surfaces. In this work, we apply ab initio molecu- 
lar dynamics[3] (AIMD), which takes into account pro- 
ton dynamics and hydrogen bond network fluctuations 
in liquid water essential to acid-base reactions in small 
molecules [S- 11 1 as well as cooperative hydroxy! hydrogen 
bonding behavior specific to oxide surfaces, 12|, |l3| to in- 
vestigate the enigmatic SiOH deprotonation equilibrium 
constant as a function of structural motifs. 

Measurements of interfacial pKa (defined as — logj^Q 
Ka, where Ka is the acid dissociation constant) have 
been revolutionized by surface-sensitive second harmonic 
generation (SHG) and sum frequency vibrational spec- 
troscopy (SFVS) techniques. [ij, [11 In 1992, Ong et 
al.fl^ demonstrated that 19% of silanol groups on fused 
silica surfaces exhibit a pKa of 4.5, about the same 
as vinegar (acetic acid), while 81% exhibit pKa=8.5. 
SFVS experiments on a-quartz reached similar conclu- 
sions and further suggested that the low- acidity silanol 



groups reside in regions with strong water-water hydro- 
gen bonds. llTl A titration study on silica gel (amor- 
phous silica) [l8| and X-ray photoelectron spectroscopy 
measurements on quartz|19| also independently demon- 
strated the existence of SiOH groups with pKa between 
4 and 5.5. Such qualitative agreement on different forms 
of silica is expected because liquid water is known to 
react with even crystalline silica to form an amorphous 
layer. [^[23| These measurements suggest that the earlier, 
single pKa ~6.8 reported in amorphous silica titration 
experiments [21j may reflect a composite of two types of 
SiOH. 

The acidities of surface silanol groups have been as- 
signed to different chemical connectivities or inter-silanol 
hydrogen bonding. In accordance with the literature, 
we differentiate SiOH groups according to whether they 
are directly hydrogen bond to other SiOH ("H-bonded," 
Fig. [T^), or are not so hydrogen-bonded ("isolated," 
Fig- lib); whether the 4-coordinated Si atom of the SiOH 
is part of 3 covalent Si-O-Si- hnkages ("Q^," Fig. [ij;), 
or only part of 2 Si-O-Si ("Q^," Fig. [TJi). It has been 
suggested that the ratio of H-bonded to isolated SiOH 
is about 1 to 4, similar to the relative occurrence of 
pKa=4.5 and 8.5;[l^ thus pKa=4.5 has been ascribed 



to isolated silanol groups. fl6l |22h24| | On the other hand, 
the Q^:Q^ ratio has also been described as either ap- 
proximately 1:4 or 4:1, which has prompted assignment 
of the pKa=4.5 SiOH group to either (Refs.[25[[2i) or 
Q"^ (Refs. [27i - [29l ). The conflicting estimates of the ratios 
of silanol groups with different structural motifs [25i - [29| 
likely reflect the difficulty brought about by the tendency 
of liquid water to react with crystalline silica. [2^ One in- 
disputable experimental flnding is that the silanol surface 
density is crsiOH ~ 4.6 nm~^ on well-soaked amorphous 
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FIG. 1: (a)-(d) Four types of SiOH groups discussed in this 
work, (a) Hydrogen-bonded; (b) isolated; (c) Q'^; (d) Q^. 
Oxygen, silicon, and hydrogen atoms are colored in yellow, 
red, and white, respectively. 



samples. 121 |3C1 On the theoretical side, static geochemical 
models. [3 ll. |32| which do not account for aqueous phase 
hydrogen bonding and dynamical proton motion, have 
also been applied, but they have not yet explained the 
two observed acidity constants, 16l4l8l| while quantum 
chemistry or DFT methods with a dielectric continuum 
treatment of the bulk water environment have been lim- 
ited to calculating the pKa of small silica fragments, [s^- 
IsBi] DFT modeling of amorphous silica slabs have also 
been considered, [3a| but the pKa estimates therein of- 
ten do not treat water explicitly or dynamically. (See 
Supporting Information (SI) for discussions of the signif- 
icance of hydrogen-bond network fluctuations and excess 
proton hopping.) 

AIMD simulations have successfull y r eproduced the 
pKa of molecules in aqueous solution[8l4lC| and should be 
particularly well-suited for distinguishing relative SiOH 
pKa that are 4 pH units apart in different environments, 
provided we can demonstrate that reproducible pKa for 
chemically equivalent SiOH's can be predicted. Coupled 
with static high-level quantum chemistry corrections, 
they provide the most rigorous predictions for liquid state 
reactions. As computing power has increased, AIMD 
modeling of liquid water- material interfaces has become 
viable. [l3, 33-4l| although it remains costly because 
water dynamics is slower at interfaces. 42, [i^ Further- 
more, investigation of multiple reaction sites and/or crys- 
talline facets is often necessary when dealing with ma- 
terial surfaces. In this work, we study the bimodal 
acid-base behavior of silanol groups [l6|-[l8j by perform- 
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FIG. 2: The heterogeneous SiOH environments examined in 
this work, (a) Hydroxylated /3-cristobalite (100) surface. The 
SiOH groups have crsiOH ~ 8 nm~^, are and H-bonded. 
pKa=7.6 ± 0.3. (b) Hydroxylated ^-cristobalite (100) sur- 
face with one SiOH group replaced with a SiH to break the 
chain of hydrogen bonds. The tagged (deep blue) SiOH (Q^ 
and isolated) exhibits pKa=8.9±0.3. (c) Reconstructed /3- 
cristobalite (100) surface, asiOH ~ 4 nm~^, Q"^ and isolated; 
pKa=8.1±0.5 (6 layers of water) and 7.0±0.4 (4 layers), (d) 
The structure in panel (c) comes from removing atoms shown 
here in blue, and attaching the resulting undercoordinated 
Si and O atoms. Panels (a)-(d) represent ~1.5 simulation 
cells in the lateral direction, (e) (H3SiO)3SiOH, which is 
and isolated, exhibits pKa=7.9±0.5. (f) Top half of a recon- 
structed quartz (0001) surface model containing cyclic silica 
trimers (Si-0)3, asiOH ~ 2.3 nm^^, Q"^ and H-bonded; one 
SiOH resides on a trimer ring (pKa=5.1±0.3), the other does 
not (pKa=3.8±0.4 and pKa=4.8±0.4 depending on whether 
a nearby trimer ring breaks; see text). Si, O, and H atoms are 
in yellow, red, and white, respectively, (b) and (e) are finite- 
temperature AIMD snapshots, with water molecules omitted 
for clarity, while (a), (c), and (f) are shown at T=0 K. 



ing AIMD simulations to directly calculate the pKa value. 
Given the absence of well-defined water-crystalline silica 
interfaces, '231 and the fact that the precise atomic struc- 
ture of amorphous silica surfaces is unknown, we examine 
six distinct, representative silanol environments. These 
include hydroxylated /3-crystobalite (100) (Fig. [5^), hy- 
droxylated /3-cristobalite (100) with one SiOH removed 
(Fig. [2)d), reconstructed /3-cristobalite (100) (Fig. [2]:), a 
molecular system (Fig.[2]3), and two distinct SiOH on re- 
constructed quartz (0001) (Fig. [2f). They represent the 
SiOH motifs proposed to be responsible for pKa=4.5 or 
8.5 in the literature (Figs. [T^-d). 



COMPUTATIONAL METHODS 

AIMD simulations apply the Perdew-Burke-Ernezhof 
(FEE) functional . 1441 the Vienna Atomic Simulation 
Package (VASP),gl El] a 400 eV energy cutoff, F- 
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point sampling of the Brillouin zone, deuterium mass 
for all protons, and a 0.375 fs time step at each Born- 
Oppenheimer dynamics time step. The trajectories are 
thermostat at T—425 K; elevated temperature is needed 
to represent liquid water properties when the PEE func- 
tional is applied and quantum nuclear effects are ne- 
glected, which is the case herein (SI, Sec. SI). Four 
to six umbrella sampling windows of 20 ps production 
trajectory length each are used per pKa calculation on 
hydroxylated /3-cristobalite (100) (Figs. [2^, b) and recon- 
structed /3-cristobalite (Fig. [St) surfaces. These periodi- 
cally replicated simulation cells measure 10.17x10.17x26 
A'^ . Their lateral dimensions are commensurate with sim- 
ulation cells often used for AIMD studies of pure water 
structure or ion hydration. This is one of the reasons we 
have focused on the (100) rather than the (111) surface 
of /3-cristobalitc, which actually has a csiOH similar to 
that of amorphous silica. [3| The more computationally 
costly reconstructed quartz (0001) system has a cell size 
of 10.0x17.32x24 A^, and 12 to 18 ps trajectories are 
used per window. For each crystalline silica simulation 
cell, we always start from crystal slab structures opti- 
mized at zero temperature using Density Functional The- 
ory (DFT) and the PBE functional. [44| Next we switch to 
the CHARMM SiOa force ficld[47] and the SPC/E model 
for water. [4^ The number of water molecules occupying 
the simulation cell is determined using these force fields, 
the Grand Canonical Monte Carlo (GCMC) technique, 
and the Towhee code.[4§| With this approach, the spaces 
between hydroxylated and reconstructed /3-cristobalite 
surfaces are filled with 58 and 63 water molecules, re- 
spectively, which amount to roughly six layers of wa- 
ter, while the reconstructed quartz simulation cell con- 
tains about 63 (about 4 layers of) water molecules. The 
(H3SiO)3SiOH pKa calculation utilizes a (12.42 A)^ cell 
with 57 H2O and 20 ps sampling trajectories. Finally, to 
check system size dependences, umbrella sampling sim- 
ulations for a smaller reconstructed /3-cristobalite (100) 
simulation cell, 20 A in the z-direction and containing 
44 (4 layers of) water molecules, are also conducted for 
at least 10 ps per window. 

pKa has been reported for molecules in liquid water 
using the AIMD technique. 8 10] It is related to the stan- 
dard state deprotonation free energy AG'") via — log^p 
exp(-/3AG'(°)), where /3 is 1//cbT and 
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dRA{R) exp[-/3AT^(i?)] \ .(1 



Here Go denotes 1.0 M concentration, R is the reaction 
coordinate, A{R) is a phase space factor to be discussed 
below, i?cut is the cutoff distance delimiting the reac- 
tion and product valleys in the free energy landscape, 
and W{R) is the potential of mean force which pro- 
vides the information needed to compute the free energy 
of deprotonation. Regardless of the reaction coordinate 



used, [1,1^ W{R) generally do not exhibit turning points 
in the deprotonated region, and -Rcut can be taken as the 
onset of the plateau where W{R) — > 0. 

The umbrella sampling method [50] is used to compute 
the W{R) associated with SiOH deprotonation. A four- 
atom reaction coordinate R (Fig. [3^) is found to work 
best under our simulation conditions. It controls what 
we call the "wandering proton" problem. We label the 
first, second and third neighbor H2O molecules of the 
SiO^ oxygen (green sphere) shown in Figs. [3}D-d "water 
1" (O depicted red), "2" (blue), and "3" (pink). When 

R 0.4 A (Fig. Eb), the SiOH bond is intact. As R 

decreases to —1.1 A (Fig. [St), the SiOH proton is trans- 
ferred to a "water 1" which has been hydrogen-bonded 
to the SiOH group, yielding a SiO~-H30+ contact ion 
pair. As R further decreases (Fig. [Sji), a proton origi- 
nally residing on "water 1" is now transferred to a second 
water molecule ("water 2"), creating a water-separated 
SiO^/HsO"^ pair, at which point the deprotonation re- 
action is almost complete. This analysis appears consis- 
tent with insights from a transition path sampling AIMD 
simulation, lllj] as follows. The Fig. configuration, 
with the excess proton and the SiO~ separated by two hy- 
drogen bonds, form a possible free energy dividing surface 
between the intact and the deprotonated acid species, 
provided that the electric polarization between the two 
states, arising from the surrounding water molecules, 11 
has sufficient time to equilibrate. Our simulation condi- 
tions allow such equilibration, and the excess proton is 
indeed observed to diffuse away if R fluctuates to regions 
significantly more negative than ^ —1.4 A. Umbrella po- 
tentials of the type (A/2)(i?— i?o)^ are used to sample the 
reaction coordinate R. Just as significantly, they ensure 
that R> —1.4 A and control the extent of proton trans- 
fer from "water 1" to all possible "water 2," so that at 
most a water-separated ion-pair is obtained. Otherwise, 
if the excess proton is several hydrogen bonds removed 
from the SiO^ (say if it spends a significant amount of 
time on "water 3" via proton transfer from "water 2"), 
it can start to diffuse ( "wander" ) through the simulation 
cell via the Grotthuss mechanism at 0(1) ps time scale 
per proton transfer. With tens of H2O molecules in the 
simulation cell and 10-20 ps trajectories, once the excess 
H+ leaves the second hydration shell of the SiO~ it does 
not return, and equilibrium sampling is not achieved. 
This "wandering" likely arises because the higher tem- 
perature and longer umbrella sampling trajectories than 
are generally used in the AIMD literature facilitate dif- 
fusion of the excess proton away from the SiO~. Other 
deprotonation reaction coordinates used in the literature 
are discussed in the SI (Sec. S2). They are found to give 
rise to wandering excess protons under our simulation 
conditions. As long as equilibrium sampling is achieved, 
the deprotonation free energy cost should not depend on 
the choice of coordinate. 



To apply Eq. [1] we use a method similar to Ref. B 
finding the most probable optimal 0(wator i)-H+ hydro- 
gen bond distance tq-h at each R, thus locally convert- 
ing W{R) to W{ro~ii); performing a spline fit to the 
resulting VF(ro-H); and integrating over t-q-h with a 
47rrQ_jj volume element, which takes the place of the 
phase space factor A{R) in Eq. [T] Equation [T] assumes 
that the entropic factors such as rotations of the reac- 
tant and products about the 0-H axis are adequately 
sampled in the AIMD trajectories; otherwise additional 
constraints and entropic factors are introduced, [sli Isij ] 
SI Sec. S2 shows that such constraints do not affect wa- 
ter autoionization free energies, partly because the varia- 
tion in ro-H needed to complete the deprotonation reac- 
tion is relatively small. Furthermore, we always reference 
predicted silanol pKa to that of water autoionizationQ 
computed using the same reaction coordinate and ele- 
vated temperature. This minimizes systematic error aris- 
ing from the simulation protocol (SI Sec. S2), and phase 
space contributions approximately cancel out. 

The metadynamics technique. [53l455| a promising and 
powerful alternative to umbrella sampling, has been ap- 
plied to calculate dissociation free energies on surfaces [3 7f 
and acid-base reactions of small molecules. This 
method, not yet implemented in VASP, can potentially 
be used for efficient comparative study of pKa on other 
material surfaces after it has been adapted to deal with 
the wandering proton problem. 

Gas-phase, high-level ab initio calculations are per- 
formed to check and correct chemical bonding ener- 
gies predicted with the PBE functional used in the 
AIMD simulations. The GaussianOS program suite is 
applied, [m The pertinent sample chemical reaction is 

Si(0H)4 + H2O • • • OH- ^ Si(0H)30" • ■ • H2O + H2O 

_ (2) 

Geometries are optimized, and harmonic vibrational fre- 
quencies computed, with density functional theory using 
the B3LYP method [13, [H^ and the 6-m++G{d,p) ba- 
sis set. At the B3LYP geometries, energies are computed 
with the coupled-cluster singles and doubles method in- 
cluding a perturbative correction for triple substitutions, 
CCSD(T),[5^ using the aug-cc-pVDZ basis set. Basis set 
incompleteness corrections are added to the CCSD(T) 
energies. Finally, zero-point vibrational energy correc- 
tions computed from the B3LYP/6-311-|--|-G((i, p) fre- 
quencies are added. Using this protocol, the high-level 
ab initio calculation yields a reaction energy of -30.51 for 
Eq. [21 Gas-phase VASP-based PBE calculations, con- 
ducted with an energy cutoff identical to that in AIMD 
simulations, predict -27.17 kcal/mol. These numbers do 
not include the 2.24 kcal/mol zero point energy (ZPE) 
corrections. Thus, the overall AIMD reaction energies 
should be corrected by a modest -1.10 kcal/mol. The 
basis set extrapolation procedure may have a systematic 
error larger than 1 pH unit (SI Sec. S3) but this does not 
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FIG. 3: (a) The 4-atom reaction coordinate R, illustrated for 
silicic acid in water but is similar for all silanol containing 
species. Panels (b)-(d) are snapshots from AIMD deproto- 
nation simulations, with outershell H2O molecules removed 
for clarity reasons. As deprotonation proceeds, R progresses 
from intact SiO-H {R ~ -0.4 A, panel (b)) to SiO" H3O+ 
contact ion pair (i? ~ —1.0 A, panel (c)) and then, via a Grot- 
thuss proton transfer, to a solvent-separated SiO'/HsO"'' pair 
{R ~ —1.32 A, panel (d)). Yellow, red, white, and green 
spheres represent Si, O, H, and the "0~" atoms, respectively. 
Additionally, the water O atoms which are second and third 
nearest neighbor to the SiO~ oxygen are colored blue and 
pink, respectively. pKa predictions for silicic acid will be re- 
ported in a future publication. 



affect the relative pKa of different SiOH groups. 

See the SI for details about constraints introduced to 
prevent proton attacks on SiO~, the AIMD intializa- 
tion protocol, further justifications for adopting the four- 
atom reaction coordinate, extrapolating quantum chem- 
istry results to the infinite basis set limit. 



RESULTS 

Hydroxylated 13-cristobalite (100): chemically homoge- 
neous SiOH. We first show that our simulation protocol 
predicts reproducible pKa for chemically equivalent SiOH 
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groups on the hydroxylated /?-cristobalite (100) surface 
(Fig.jJk)- This wen-studied model crystaUine surface ex- 
hibits crsiOH=8 nm"^, larger than the experimental value 
of 4.6 nm^^ for amorphous sihca.pj At zero-temperature, 
it features two types of silanol groups: alternating hy- 
drogen bond donors and acceptors arranged in chains 605 
(Fig. [2^) not found in smah molecules. 8-10] This fea- 
ture is dynamically preserved in our finite-temperature, 
aqueous-phase simulations (SI Sec. 84). 

Figure |4^ shows that two chemically equivalent, hy- 
drogen bond-accepting SiOH groups on this surface 
are predicted to exhibit deprotonation W{R) within 
0.5 kcal/mol of each other. At large negative values of 
the reaction coordinate R, the deprotonated SiO~ is sta- 
bilized with three hydrogen bonds (i.e., iVu,=3, Figs. |4)d 
anddt). At R -0.8 A, the SiO-H bond is only par- 
tially broken, iV^, < 3, and the local W{R) is not sensi- 
tive to the slight difference in iVy, in the two simulations 
that arises from statistical noise. These observations ap- 
pear consistent with a recent two-dimensional potential 
of mean force analysis of formic acid deprotonation. [To] 
Accounting for zero-point energy, correcting the AIMD 
functional with more accurate quantum chemistry meth- 
ods, and referencing Eq. [I]to the water autodissociation 
constant pKw=14,^j we estimate pKa values of 7.5 and 
7.7, close to the less acidic pKa value reported by Ong et 
aL[16j The standard deviation is 0.3 pH unit. Multiple 
deprotonation on this surface is discussed in SI Sec. S5. 

Heterogeneity: Isolated, H-bonded, Q^, and SiOH 
all exhibit pK^ > 7. We next show that, contrary to pre- 
vious hypotheses, 16, 23-2^ isolated, H-bonded, Q^, and 
Q'^ silanol groups all exhibit pKa > 7.0. We first create 
an isolated silanol group by replacing a hydrogen bond- 
donating SiOH group on the hydroxylated /3-cristobalite 
(100) surface with a SiH so that its neighboring SiOH 
group is no longer H-bonded (Fig. [5}d). Figure [S] shows 
that this isolated SiOH exhibits pKa=8.9±0.3, and is less 
acidic by 1.2-1.4 pH unit than when the SiOH hydro- 
gen donor is present (Fig. [2^). This is entropically rea- 
sonable because a hydrogen bond-donating SiOH part- 
ner stabilizes the neighboring SiO^ alongside two water 
molecules, while three water molecules are required for 
an isolated SiO^. AIMD correctly accounts for this effect 
because it models H2O and SiOH on the same dynami- 
cal footing and because water- water and water-silanol hy- 
drogen bond energies are similar. [601] Hydroxy Is on oxides 
with more ionic character than Si02 form stronger hydro- 
gen bonds, and indeed the relative abundance of inter- 
hydroxyl hydrogen bonding may partially be responsible 
for the crystal facet-dependent acidity of Q;-Al9 0-i.[6ll| 
This will be the subject of future, comparative studies. 

To compare and Q'^ silanol groups, we reconstruct 
the hydroxylated /3-cristobalite (100) surface by condens- 
ing every other pair of H-bonded SiOH groups into a 
SiOH and a H2O molecule (Fig. [2t). This involves 



the elimination of a hydrogen-bond donating OH group 
plus the proton on its adjacent, hydrogen-bond accepting 
SiOH (Fig.[2Ji). The resulting undercoordinated Si and O 
atoms are joined together to form a covalent bond, in the 
process pulling apart the remaining H-bonded SiOH pairs 
so they are now isolated from one another. A similar 
structural motif has been considered in the literature. 
This surface has crsiOH=4 nm^^, with all SiOH groups 
being Q'^, isolated, and residing on silica rings contain- 
ing at least 5 Si atoms. Such rings should be unstrained, 
unlike 3-member (Si-0)3 rings discussed below. The pKa 
is found to be 8.1±0.5 (Fig.©. 

We also consider a (SiH3)3SiOH molecule featuring 
an isolated, Q'^ SiOH group (Fig. [2^), which exhibits 
a comparable pKa=7.9±0.5. Thus, in general, Q'^ and 
silanol groups do not exhibit pKa's that differ by 4 
pH units as previously proposed. |23l - l29| Over the range 
4 nm~^ < (TsiOH < 8 nm~^, the precise value of csiOH 
has little effect on pKa. 

To some extent, all our crystalline silica models are 
nano-slits with thin water slabs confined between the 
surfaces. As the water content decreases, the dielectric 
solvation of SiO~ and H+ species should decrease, while 
intact SiOH groups should be weakly affected. Thus one 
expects a lower acidity and a higher pKa in strongly 
nano-confined aqueous media. [63j| To examine confine- 
ment effects, a pKa calculation is performed for a smaller 
reconstructed /3-cristobalite (100) simulation cell, 20 A in 
the z-direction, containing 4 layers of water. This system 
actually yields a lower pKa=7.0±0.4 — but is lower only 
by 1.1 pH units (dashed brown line in Fig. [5]). It is pos- 
sible the unexpected pKa decrease between the 6- and 4- 
water layer models arises from anomalies in the hydrogen 
bonding network not apparent from visual inspection of 
water configurations. The water density is also affected 
by the confinement. [64] In the 4-layer model, all water 
molecules are at most two layers away from the crys- 
talline silica surfaces, and AIMD conducted with GCMC- 
predicted water content is found to yield a second layer 
H2O density that is 18% above 1.0 g/cc. However, this 
is unlikely to change the dielectric response sufficiently 
to lower the pKa by 1 unit. Assuming one can apply the 
Born hydration formula for excess proton hydration in 
this heterogeneous medium, AGhyd ^ X(l — l/eo) where 
X PS —264 kcal/mol for the proton. [65] If the 6-layers 
of water have eo=80, the 4-layer system must exhibit 
£0=140 to make hydration more favorable by 1 pH unit 
when in fact confinement generally reduces the dielec- 
tric constant of water. [l2| In any case, the discrepancy 
in pKa is actually within two standard deviations and 
may simply arise from statistical uncertainties. This test 
suggests that confinement effects are not large for the slit 
pore geometry^] down to about 1 nm slit widths. There- 
fore our reported pKa for 6-water layer systems should 
be good approximations of pristine crystalline silica sur- 
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FIG. 4: Deprotonation free energy on hydroxylated /3- 
cristobalite (100) (Fig. [2^). (a) Potential of mean force 
W{R). The most negative R values represent SiO~ while 
R >r^ —0.8 A indicates intact SiO-H bonds. Red and black 
lines: two different but chemically equivalent SiOH groups; 
blue: water autoionization. (b) Mean hydration number, N^, 
in each umbrella sampling window, defined as the number 
of protons within 2.5 A (a typical hydrogen bond distance) 
of the SiO~ oxygen, (c) Snapshot of the interface between 
water and hydroxylated /3-cristobalite (100), replicated three 
times in the lateral direction, (d) Side view of slab model 
with silanol groups forming hydrogen bond chains, highlight- 
ing a Q^, H-bonded SiO~ group in its representative hydrogen 
bonding environment. Yellow, red, and white spheres depict 
Si, O, and H atoms, respectively. 



faces in contact with bulk liquid water. We also stress 
that all slab geometries in Figs. [5^-c have been studied 
using 6-water layer simulation cells, and their relative 
pKa should be mostly free of system size effects. How- 
ever, two-dimensional confinement in cylindrical amor- 
phous silica nanopores [12n23i . i66H68j may have a stronger 
impact on pKa.|63j 

High acidity and chemical reactions on strained recon- 
structed quartz surfaces. Finally, it seems imperative to 
demonstrate the possibility of an unusually acidic silanol 
group. The following "computational existence proof" is 
necessarily more speculative than the conclusions about 
Q^, Q^, isolated, and H-bonded SiOH pKa discussed 
above, but it emphasizes the likely role of defected re- 
gions when accounting for pKa=4.5. 



Having considered crsiOH=8 and 4 nm^^, we examine 
an even lower SiOH surface density. A recent experi- 
mental study has attached crystal violet dyes to deproto- 
nated silanol groups. Based on the flat, 120 nm^ surface 
area of the dye molecule, it is proposed that strong acid- 
ity is correlated with local SiOH surface density csiOH 
< 0.83 nm~^.^22i] In our DFT calculations, the opti- 
mal geometry of crystal violet when covalently bonded 
to (HO)3SiO~ is not flat but is substantially distorted. 
However, this suggestion of low SiOH surface density be- 
ing associated with the more acidic SiOH appears consis- 
tent with other experimental and theoretical observations 
discussed below. 

Since most cuts through crystalline forms of silica yield 
surfaces with substantial hydroxylation, [69] we investi- 
gate a reconstructed, completely dehydroxylated quartz 
(0001) model, featuring (Si-0)3 trimer rings, predicted 
to be metastable in vacuum. [toI [Tlj 2.3 SiOH groups per 
square-nanometer are re- introduced by removing two sur- 
face Si atoms and performing further reconstruction and 
hydroxylation to keep all atoms fully coordinated. This 
yields two types of silanol groups which are hydrogen- 
bonded to each other; one member of the pair resides 
on a cyclic trimer while the other does not (Fig. [2f). 
Regions devoid of SiOH and dominated by siloxane (Si- 
0-Si) bridges are hydrophobic, and this model surface 
may therefore be consistent with hydroxyl groups in hy- 
drophobic pores reported to be unusually acidic on other 
oxides, [t^ We conduct one deprotonation umbrella sam- 
pling simulation of a silanol group residing on a cyclic 
trimer (Fig. [SJ solid violet curve), and two simulations 
where the SiOH docs not reside on a trimer (dashed 
and dot-dashed violet). Figure [5] indeed shows that 
these SiOH groups exhibit pKa 5.1±0.3, 4.8±0.4, and 
3.8±0.4, respectively — close to the experimental value of 
4.5. More significantly, their average pKa are separated 
from the median of all other SiOH groups previously ex- 
amined in this work by 3.4 pH units. 

Unlike cyclic silica tetramers or larger Si-0 rings, cyclic 
silica trimers are strained. [lIMIzl At zero temperature, 
in the absence of water, the Si atom of the SiOH group 
residing on a silica trimer ring exhibit Si-O-Si angles of 
137.5°, 130.4°, and 132.0°. For the SiOH group not re- 
siding on a silica trimer, the angles are 154.9°, 131.8°, 
and 147.8°. A few of these angles deviate substantially 
from the ideal, unstrained Si-O-Si value of approximately 
145°. This likely accounts for the low pKa computed for 
these SiOH groups. Other trimer rings not decorated 
with SiOH are also strained, and the slow spatial decay 
of their surface strain fields 73| may also contribute to 
the high acidity of SiOH group not residing on them. 

Within hours in moist air,[l[ [tI, [t^ cyclic trimer- 
containing surfaces are known to incorporate water and 
break open to reduce strain and increase the local asiOH- 
At the water-reconstructed quartz interface, during um- 
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R(A) 

FIG. 5: Distribution of pKa. Black and red lines: W{R) 
for SiOH on the hydroxylated /?-cristobalite (100) surface 
(Fig. [2^); green: same surface but with hydrogen bond- 
donating SiOH artificially removed (Fig. [Us); brown, solid 
and dashed: reconstructed /3-cristobalite (100) (Fig[2j;), with 
6 and 4 layers of water in the simulation cell, respectively; 
yellow: (H3SiO)3SiOH (Fig. [2^); violet: reconstructed quartz 
(0001), (Fig.O), with the SiOH residing on a silica trimer ring 
(solid) or otherwise (dashed and dot-dashed). Blue: W^(-R) for 
water autoionization. 



brella sampling deprotonation of the SiOH group residing 
on a 3-member ring (Fig. [2f), we indeed observe a water 
molecule forming a transient bond with another Si atom 
on a Si-0 trimer 6 A away from the tagged SiO~ within 
picoseconds (Fig.jBb). The resulting 5-coordinated Si has 
been observed in simulations (tH [tgI. \71\ and found to be 
the intermediate in the trimer ring-opening mechanism 
on wet silica amorphous surfaces in reactive force field 

[liS] Our AIMD 



and molecular orbital calculations, 
trajectories show that this mechanism remains operative 
at explicit liquid water-silica interfaces; proton hopping 
via the Grotthuss mechanism occurs readily, enabling 
the H2O adsorbed on the surface Si to lose a proton to 
bulk water, forming a new SiOH group within ^ 10 ps 
(Fig. int). Then, in 2 of the 4 sampling windows, a Si- 
O bond on the now 5-coordinated Si breaks to open the 
OH-incorporated (Si-0)3 ring and irreversibly introduce 
another new SiOH (Fig. |6ji). Our trajectories thus dif- 
fer from a recent molecular dynamics study of the liquid 
water-reconstructed quartz interface, where the adsorbed 
water molecule, not described by a reactive force field, ul- 
timately desorbs from the 5-coordinated surface Si atom 
without inducing chemical reactions. 71 



These irreversible side reactions prevent strict equi- 
librium sampling needed for W{R) calculations. Fortu- 
nately, for the SiOH residing on a cyclic trimer (Fig. [S]) , 
analysis of the pre- and post ring-breaking statistics re- 
veals that the nearby chemical reaction has little effect on 
its pKa. We further analyze trimer ring-opening effects 
on the pKa of another SiOH group, this one not residing 
on a surface 3-member ring. A H2O incorporation reac- 
tion also occurs in the neighborhood of this tagged SiOH 
(Fig.[7]). We split the sampling windows into two groups: 
(A) those without irreversible hydrolysis of a nearby sil- 
ica ring (Fig.[7K); and (B) those with trimer ring breaking 
and formation of two new SiOH groups (Fig. Uh)- Then 
two complete sets of sampling windows spanning the en- 
tire deprotonation pathway are spawned from these seed 
windows, yielding two pKa: case A, pKa=3.8 (Fig. [5j 
dashed violet curve); and case B, pKa=4.8 (Fig. [Sj dot- 
dashed violet). The results show that H2O incorporation 
and a single ring-breaking event nearby does appear in- 
crease the pKa of the tagged SiOH not residing on a 
trimer ring. However, the increase is only 1.0 pH unit, 
almost within statistical uncertainties. In case (A), the 
three Si-O-Si angles on the silica trimer ring average to 
121.1°, 134.0°, and 132.2° along the trajectory; the first 
refers to the angle where both Si are below the silica 
surface ("buried"). In case (B), this angle linking the 
buried Si atoms relaxes significantly to 140.5°. The sec- 
ond angle, which involves the surface Si and a buried Si, 
remains almost unchanged at 136.0°. (The third linkage 
is destroyed during ring opening.) Si atoms which no 
longer participate in strained Si-O-Si linkages should be 
more stable against H2O attack. 

From these considerations, we conclude that, despite 
interference from H2O incorporation reactions, there re- 
mains a statistically significant difference in the pKa's 
on this surface, and the pKa's ranging from 7.0 to 8.9 
for all other silanol groups investigated before (Fig. [5]) . 
This finding suggests that strain, low local silanol surface 
density, hydrophobicity, and low pK^ are correlated on 
amorphous silica surfaces. Indeed, atomistic model sur- 
faces with low local <7R inH regions almost always exhibit 
3-member rings. 7^ 76, 7^ These regions may be in dy- 
namic equilibrium with solvated silica fragments in solu- 
tion, constantly being dissolved/hydrolyzed and reconsti- 
tuted when dissolved fragments re-nucleate on hydroxy- 
lated regions. [Tli Isol - fs^ The dynamic equilibrium has re- 
cently been demonstrated in Monte Carlo simulations [8lj| 
using a reactive silica force field, js^ 

We have only observed one ring-opening reaction on 
each surface. The limited AIMD trajectory length does 
not conclusively allow us to predict how many trimer 
rings persist at the liquid water-amorphous silica inter- 
face as a function of time. Therefore we do not defini- 
tively assign this structure to the observed pKa=4.5 
SiOH group, and instead pose it as a challenge to experi- 
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FIG. 6: Water incorporation onto reconstructed quartz. 
Initial equilibrated snapshot of SiO~ contact ion pan 
on the reconstructed quartz (0001) surface where the SiO" 
resides on a S-member ring (truncated in this figure). The 
tagged SiO~ oxygen is in deep blue; its H-bond SiOH donor 
is out of the frame. All other Si, O, and H atoms are in yel- 
low, red, and white, respectively, (b) Snapshot after ~ 3 ps. 
To the left of the SiO~, a water molecule (light blue) has at- 
tached itself to a Si atom on another silica trimer ring, which 
is unhydroxylated. That Si becomes 5-coordinated. (c) This 
water molecule loses a proton to bulk water; one of the Si- 
O bonds not on the trimer ring becomes stretched. This is 
reversible as long as another SiOH is not created, (d) In an- 
other sampling window for this SiOH, the trimer ring breaks 
open instead, forming a new SiOH after extracting a proton 
from water. 



mental work, including single molecule spectroscopy, 84 
to determine whether they are sufficiently abundant over 
time to account for the 19% of all silanol groups shown 
to exhibit high acidity. [3, [13 We also point out that, 
while cyclic silica trimers are well known to react with 
moist air, other popular crystalline silica model surfaces 
should also be hydrolytically unstable, [i^l Thus, the hy- 
droxylated a-quartz (0001) and /3-cristobalite (100) sur- 
faces, with CTsiOH ^ 9 and 8 nm~^ respectively, are of- 
ten used as models to study the interface between liquid 
water and generic silica solids in classical MD simula- 
tions. However, if these simulations permit chemical re- 
actions over long enough times, we speculate that some 
of the SiOH groups on such surfaces may also react with 
water 2^ 81 1 in a way to reduce the silanol surface density 
towards the amor pho us silica crsiOH==4.6 nm~^ observed 
in experiments. 0,|33| 

Discussions AIMD-based potential of mean force cal- 
culations simulations have been demonstrated to yield re- 
producible pKa for chemically equivalent silanol groups. 




FIG. 7: Trimer silica ring and Si-O-Si angles. Configura- 
tions are taken from AIMD pKa calculations where the SiOH 
group involved in the pKa calculation does not reside on a 
silica trimer ring. That tagged SiOH is off- frame; all H2O 
molecules are removed for clarity, (a) Snapshot along AIMD 
trajectory where the 3- member ring remains intact (case A). 
The 3 Si-O-Si angles within the ring average to 121.1° (be- 
tween Si atoms on the trimer ring colored yellow and dark 
blue), 134.0° (yellow/pink), and 132.2° (pink/dark blue), (b) 
Incorporation of a H2O (its oxygen colored light blue) opens 
the ring (case B). The surviving Si-O-Si angles are 140.5° 
(yellow/dark blue) and 136.0° (yellow/pink). O and H atoms 
are represented by red and white spheres respectively. 



The statistical uncertainties of our simulation protocol 
are estimated to be about 0.3-0.5 pH unit, consistent 
with explicit calculations on two chemically equivalent 
SiOH. Therefore these simulations should reliably distin- 
guish relative pKa of heterogeneous SiOH groups 4 pH 
units apart. Resolving hydroxyl pKa on other surfaces 
may remain a challenge if the acidities are less widely 
separated. 

Our pKa calculations suggest that comparative stud- 
ies between Si02 and other material surfaces will be ex- 
tremely interesting. One intriguing candidate surface is 
a quartz surface densely functionalized with carboxylic 
acid groups. To our knowledge, this is the only other 
material surface which clearly exhibits bimodal pKa 
behavior, [s^ Other candidates are the different facets of 
crystalline alumina, which may feature several pKa 's un- 



resolved into distinct, measurable components. [72ll86ll87l | 
Structural motifs such as chemical connectivity and inter- 
hydroxyl hydrogen bonding have also been invoked to ex- 
plain the pKa in these systems; as mentioned in the text, 
inter-hydroxyl hydrogen bonding may affect the pKa of 
the more ionic AI2O3 surfaces more strongly than on any 
form of silica surfaces. See SI, Sec. S6, for more details 
on these material systems. Finally, SFVG spectra[lj 
and time-dependent acid-base phenomena on quartz 
can also be investigated in the future. 



CONCLUSIONS 

In this paper, we have performed AIMD pKa calcula- 
tions on five representative crystalline silica surfaces plus 
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a molecular system exhibiting different silanol (SiOH) 
structural motifs. From the results, we have conclusively 
shown that the more acidic of the two pKa observed in 
experiments cannot, as previously proposed. (22n2i^] be 
explained by the existence of silanol groups with certain 
chemical connectivities or inter-silanol hydrogen bond- 
ing. In fact, we find pKa ^ 4.5 silanol groups only on 
strained surfaces with sparse silanol coverage. While our 
demonstration of the existience of such low pKa SiOH 
groups is necessarily somewhat speculative, this study 
highlights the role of defected regions as the most promis- 
ing candidate to explain the elusive bimodal acid-base 
behavior of silica surfaces. Assigning structural 

motifs to the more strongly acidic SiOH groups is par- 
ticularly crucial in non-reactive force field-based mod- 
eling of silica nanofiuidic channels, where the preferen- 
tially deprotonated SiOH sites at neutral pH have to be 
assigned in a static way. 12, 23, 66, 67, [s^ In the 
process of studying the acid-base behavior, we also ob- 
serve irreversible, water-assisted ring-opening reactions 
of strained silica trimer rings in contact with liquid wa- 
ter, The reaction was previously studied on wet silica 
surfaces: [76. 76 . 7^ our AIMD simulations demonstrate 
that a similar mechanism is operative at liquid water- 
silica interfaces. 
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